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Abstract 

This paper introduces a new method for proving global stability of fluid flows through 
the construction of Lyapunov functionals. For finite dimensional approximations of fluid 
systems, we show how one can exploit recently developed optimization methods based 
on sum-of-squares decomposition to construct a polynomial Lyapunov function. We then 
show how these methods can be extended to infinite dimensional Navier-Stokes systems 
using robust optimization techniques. Crucially, this extension requires only the solution 
of infinite-dimensional linear eigenvalue problems and finite-dimensional sum-of-squares 
optimization problems. 

We further show that subject to minor technical constraints, a general polynomial Lya- 
punov function is always guaranteed to provide better results than the classical energy 
methods in determining a lower-bound on the maximum Reynolds number for which a 
flow is globally stable, if the flow does remain globally stable for Reynolds numbers at least 
slightly beyond the energy stability limit. Such polynomial functions can be searched for 
efficiently using the SOS technique we propose. 
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1. Background and problem statement 



In this paper we propose a new analytical method for determining whether a fluid flow is 
globally stable. This new approach has its origins in two hitherto distinct research areas. 
The first of these is the classical energy approach of [H [16] , which provides conservative 
lower bounds on the stability limits of flows by analyzing the time evolution of the energy 
of flow perturbations. The other is the emerging field of sum-of-squares (SOS) optimization 
over polynomials, which can be used to prove global stability of finite-dimensional systems 
of ordinary differential equations with polynomial right-hand sides [131 [H] ■ It is our hope 
that the present text is written in such a way that it will be understandable to researchers 
from either of these two areas, which until now have remained almost completely isolated. 
We are aware of only two other publications where SOS methods have been used to analyze 
the behavior of similar systems governed by partial differential equations, namely |T1, 20J. 

We apply SOS methods to the problem of assessing global stability of an incompressible 
flow. The velocity w and pressure p' of a flow of viscous incompressible fluid, evolving inside 
a bounded domain f2 with boundary d£l under the action of body force f , is governed by 
the Navier-Stokes and continuity equations 

+ w • Vw = -Vp' + — V 2 w + f (la) 
at Re 

V-w = 0, (lb) 

with a boundary condition w = on d£l. Here Re is the Reynolds number, which is a 
dimensionless parameter indicating the relative influence of viscous and inertial forces in 
the flow. In what follows we will make extensive use of an inner product of vector fields 
defined as 

(u, v) := / u • v dfl 
Jn 

with the usual £2 norm |.| defined as ||u|| 2 := (u, u). Similarly, we define (using standard 
Einstein summation notation throughout) 



(Vu ' Vv> S -J a ww ia 



and || Vu]| := (Vu, Vu). 

We say that a steady solution w = u, p' = p of the system ([I]) is globally stable if, for each 
e > 0, there exists some S > such that ||w — u|| < 5 at time to implies that ||w — u|| < e 
for all time t > to. We say that it is globally asymptotically stable if in addition w — > u 
as time t — > 00 for any initial conditions. For the system ([!]), these stability conditions 
ensure laminar flow. Our principal aim is to identify the largest value Re for which these 
conditions can be guaranteed to hold for the system ([!]) . 
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The results described in this paper can be extended to other types of boundary conditions, 
most notably to the frequently encountered case of periodic boundary conditions. A useful 
property of systems with such boundary conditions is that, for any solenoidal vector field 
or fields satisfying these boundary conditions and for any scalar function <f>, the following 
hold true: (v, V0) = 0, (vi, V 2 v 2 ) = (v 2 , V 2 vi) = - (Vvi, Vv 2 ) and 

(vi, v 2 • Vv 3 ) = - (v 3 , v 2 • Vvi) , (2) 

hence 

(v,v-Vv)=0. (3) 

These properties can be proved by applying standard identities from vector calculus or sim- 
ply by integrating by parts, using incompressibility (V • v = 0) and applying the boundary 
conditions. These properties are, of course, well known. 

In particular, the identity ^ plays an important role in the energy approach to proving 
global stability [IB]. Defining velocity perturbations u := w — u and pressure perturbations 
p := p' — p, the system ([!]) can be written as 

c)w 1 

— + u • Vu + S(u, u) = -Vp + — V 2 u (4a) 
ot He 

V • u = 0, (4b) 

where <S(u, v) := u • Vv + v • Vu is introduced for compactness of notation. The time rate 
of change of the velocity perturbation energy can be obtained by taking the inner product 



of both sides of (4a) with u, to obtain the energy equation 



^£^ = ^<u,V 2 u>-(u,5(u,u)). (5) 



Note that the nonlinear term u • Vu in ( 4a ) does not feature in the energy equation because 
of the identity ^ . This is particularly useful because it allows one to obtain immediately 
an (albeit conservative) method for checking stability of the system Q, which we now 
describe briefly. 

There exists a real constant k such that for all solenoidal u, 

-^(u,V 2 u)-<u,S(u,u)) <k||u|| 2 . (6) 

We can follow the procedure in [31 p. 33-34] to find the smallest such k via solution of 
an eigenvalue problem. The smallest k satisfying ^ is the solution to the optimization 
problem 



sup 



1 "Vu|| 2 -/ n u-(Vu)-udV 



Re 



I 1 1 2 
U 
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subject to the incompressibility condition V • u = and the boundary conditions. Since 
the objective function in this optimization problem is homogeneous in u, one is free to 
optimize over the numerator only, with additional constraint ||u|| = 1. This leads to the 
eigenvalue problem 

Au=(e-^VV + Vp (?) 
V • u = 0, u\ m = 0, 

where p is the Lagrange multiplier for the incompressibility condition, A is the Lagrange 
multiplier for the unit norm condition ||u|| = 1 and e is the base- flow rate of strain tensor 
with components 

eij(x) :— 
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du l du ] 
dxi dx l 



Note the identities used in arriving at ([7]); (u, S(u, u)) = (^J^^vl) °+ (u,u- Vu) and 
u • (Vu) • u = u • e • u. Since the operator in ([7]) is self-adjoint \3\, all the eigenvalues of 
([7]) are real. If these eigenvalues are ordered by decreasing value with Ai being the largest, 
then the inequality ^ is tight with k = Ai. If the largest eigenvalue Ai < 0, then ([5|is 
always negative for ||u|| ^ and hence the energy ||u|| 2 /2 is a Lyapunov functional forMj), 
thus proving the global stability. 

A particularly nice feature of the energy approach is that proving global stability by this 
method requires only the solution of linear eigenvalue problems, even though Q is a 
system of nonlinear partial differential equations. This is a direct consequence of the 
unique advantages of using ||u|| as a Lyapunov functional, in particular the opportunity to 
exploit ([3]). In general, any other choice of Lyapunov functional would result in a stability 
condition featuring the nonlinear inviscid term u • Vu, in contrast to the more benign 
energy stability condition ([5]). 

On the other hand, the energy approach can give very conservative results, in the sense 
that the largest Re for which global stability can be proven by this method is generally 
well below the maximum Re for which the flow is generally observed to be globally stable, 
either numerically or experimentally. 

The approach proposed in the present study aims to improve this bound, using a partial 
Galerkin decomposition of the infinite dimensional system Q. Finite dimensional meth- 
ods, based on recently developed techniques in polynomial optimization, are used to define 
a Lyapunov functional that is nonlinear in a finite number of terms, while otherwise main- 
taining some of the attractive numerical advantages of energy methods for the remaining 
(infinite dimensional) dynamics. We stress that our results suggest a way of computing 
a Lyapunov functional verifying stability of the infinite dimensional system Q, and not 
some truncated finite dimensional approximation thereof. 
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1.1. Finite Dimensional Systems and the Sum- of- Squares Decomposition 

We first comment briefly on the state of the art in direct methods for computing Lyapunov 
functions for finite dimensional nonlinear systems. Suppose that the evolution of a finite- 
dimensional system with state vector a G M n is governed by a set of ordinary differential 
equations (ODEs) 

a = /(a) (8) 

with equilibrium point a = 0. We will use V a throughout to indicate the gradient of a 
scalar function defined on this n-dimensional state space, and otherwise use V to indicate 
the gradient or divergence of functions in physical space, as in ([I]). 

The origin of the system ^ is globally asymptotically stable if there exists a continuously 
differentiable Lyapunov function V : W 1 — > R such that V(0) = 0, V(a) > for all 
ael n \ {0} and V(a) = V a V ■ f{a) < for all a 6 M n \ {0} [6j. Given a Lyapunov 
candidate function V(a) and associated V(a), these conditions amount to checking global 
positivity or negativity of functions. There is no general method for performing such a 
check, nor any systematic way of constructing Lyapunov functions for general systems of 
ODEs. 

A truncated Galerkin approximation reduces the Navier-Stokes equations Q to a system of 
ODEs in exactly the form Q, but with polynomial (in fact quadratic) right hand side /(a). 
In this particular case, checking that a polynomial function V{a) serves as a Lyapunov 
function reduces to verifying the positive-definiteness of the two related polynomials V{a) 
and —V{a). However, verifying positive-definiteness of a general multivariate polynomial 
is still NP-hard in general, and is a classical problem in algebraic geometry. 

Nevertheless, there has been significant recent progress in stability analysis of polyno- 
mial systems using sum- of- squares optimization methods, which were first employed in 
the context of dynamical systems in [13] . These methods are based on a recognition that 
a sufficient condition for a polynomial function to be positive-definite is that it can be 
rewritten as a sum-of-squares (SOS) of lower order polynomial function^] Verifying this 
stronger condition, and solving other problems related to such representations, is signifi- 
cantly simpler than verifying global positivity in general. Therefore the general approach 
of sum-of-squares optimization in control applications is to search for a Lyapunov function 
V(a) and associated function V(a) that satisfy sum-of-squares conditions. 



1 This condition is not a necessary one however, apart from certain exceptional cases involving relatively 
few variables or low-order polynomials. An example of a positive-definite polynomial function that is not 
a sum-of-squares is the Motzkin polynomial y 4 z 2 + y 2 z 4 — 3y 2 z 2 + 1. A nice account of the history of this 
problem, the 17th of Hilbert's 23 famous problems posed at the turn of the 20th century, can be found 
in [E]. 
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Every polynomial V(a) of order 2k can be represented as a quadratic form of monomials of 
order less than or equal to k, i.e. in the form V = Ciji7ii(a)mj (a) . The monomials rrii(a) in 
this factorization are expressions of the form a^ 1 a^ 2 . . . a,^ n , with integer exponents ki > 
satisfying ^™ ki < k. If the matrix Cjj is positive-definite and symmetric (the latter is 
always possible), then it can be diagonalized by a suitable linear transformation of the 
monomial set. Since all the diagonal elements in the resulting expression will be positive, 
this gives a representation of the polynomial as a sum of squares of polynomials of lower 
order. Such a representation is known as a sum-of-squares decomposition. 

Hence, the problem of finding a Lyapunov function is reduced to finding coefficients = 
Cji such that this matrix is positive-definite and the corresponding matrix factorization 
representing —V(a) is also positive-definite. The relationship between the coefficients of 
V{a) and the matrix Cjj amounts to a set of linear equality constraints, with similar linear 
equality constraints relating the coefficients of V{a) and its factorization. A further set of 
equality constraints couple the coefficients of the polynomial functions V(a) and V(a). 

Problems such as that described above can be solved efficiently since the set of positive- 
definite matrices is convex. The general field of optimization theory and numerical methods 
related to such problems is known as semidefinite programming and such problems are 
solvable in an amount of time that is polynomial in the size of their problem data [21 fT8| [T9] . 
Standard software tools are freely available for posing and solving sum-of-squares problems 
[3 El H2] as semidefinite programs. 



1.2. Application of SOS methods to fluid systems 

With respect to the SOS approach, ODE systems obtained via finite-dimensional approx- 
imation of the Navier-Stokes equations require special treatment for two reasons. First, 
since the ODEs describing the dynamics are quadratic, if V(a) is of even degree (as it 
must to be positive-definite), then —V(a) is formally of odd degree, and hence will not be 
positive-definite in general. 

The second reason is more subtle. Consider the behaviour of the Lyapunov functional for 



very large values of ||u||. In this case the term S(u, u) and the viscous term 4-V 2 u in (4a) 



become small relative to the nonlinear term u-Vu, and the dynamics become approximately 
inviscid. In the limit the time derivative of a Lyapunov functional can remain negative or 
can tend to zero. In the first case the high-||u|| asymptotics of the Lyapunov functional 
would be a Lyapunov functional of the zero solution for the inviscid flow with zero forcing. 



2 Strictly speaking, the problem described here has positive-definite matrix constraints, rather than 
semidefinite constraints as in standard semidefinite programming. The conditions described above can be 
recast as semidefinite constraints via inclusion of appropriate terms, e.g. V(a) > if V(a) — e||a|| 2 = 
Cijm,i(a)mj(a) > for some small e > 0. In this case the semidefiniteness condition CV, >r is sufficient. 
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This is impossible because in such an inviscid flow the energy is conserved, i.e. the flow 
does not decay to rest. In the second case the asymptotics will be a functional remaining 
constant on any solution, i.e. it will be an invariant of the inviscid flow. 

In the 2D case there are an infinite number of such invariants, so that no conclusions can 
be drawn, but in the 3D case the only invariant known so far is energy. Hence, it is highly 
likely that the high-||u|| asymptotics of the Lyapunov functional of the viscous flow match 
||u|| itself or a monotone function thereof. Therefore, it is reasonable to limit the search for 
a Lyapunov functional to functionals with such asymptotics. Moreover, since ||w|| has the 
same high-||u|| asymptotics as ||u|| itself and since, unlike ||u|| , ||w|| decays monotonously 
in a viscous flow for any Re if ||u|| is large enough, it is sensible to seek a Lyapunov 
functional that behaves like ||w|| at large ||u||. Although this argument is not rigorous, it 
helps in guessing the structure of an appropriate Lyapunov functional. A more technical 



analysis of the asymptotic behaviour of Lyapunov functionals is given in Section 3.1 for 
the finite-dimensional case. 

We describe in Sections [2] and [3] how one can apply SOS methods to truncated ODE 
approximations to Q to obtain numerical estimates of the maximum value Re for which 
the system is global stable. We supply a numerical example illustrating the application of 
these methods in Section [6j where a stability limit approximately seven times larger than 
the value demonstrable via energy methods is obtained, and which is close to the global 
stability limit estimated by direct numerical simulation. A preliminary version of these 
results was also presented in [I]. 

There remains the question of convergence of global stability results obtained in this way 
as the number of modes in the truncated Galerkin approximation tends to infinity; global 
stability of a truncated approximation does not imply global stability of the Navier-Stokes 
solution. This problem is particularly acute since one cannot realistically expect to apply 
SOS methods to high resolution approximations of the Navier-Stokes equations, since the 
size of the related optimization problems quickly becomes unmanageable. 

We therefore demonstrate in Sections [4] and [5] how these difficulties can be overcome by 
searching for a Lyapunov functional of the Navier-Stokes system in the form V = V(a, q 2 ), 
where a is a finite-dimensional vector of the amplitudes of several Galerkin modes, and q 2 
is the collective energy of all of the remaining modes. Such an approach requires estimates 
via q for the terms stemming from the nonlinearity of the Navier-Stokes system, which is 
immediately reminiscent of the difficulties preventing proof of existence of solutions to the 
Navier-Stokes equations. In the present context it turns out, however, that the required 
estimates are available since they are needed only for the effect of higher-order modes on the 
finite set of Galerkin modes. As a result, the required estimates can be obtained by solving 
only linear eigenvalue problems in infinite dimensions and certain maximization problems 
in finite dimensions, and the resulting system can be treated using the SOS approach. 
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We show in Section |5.1| that, with a suitable basis for the Galerkin approximation, the 
proposed approach is guaranteed always to give results at least as good as the standard 
energy approach. We further show in Appendix B that if the flow remains globally stable 
in some range of Re beyond the maximum Re for which stability can be proved using the 
energy approach, then a polynomial Lyapunov function is still guaranteed to exist in at 
least part of this extended range. 



2. Finite Dimensional Flow Models 



We assume throughout that the perturbation velocity u can be written as 

u(x,t) = ai(t)ui(x) + u s (x,t), i = l,...,n (9) 

where the basis functions Uj are mutually orthogonal, solenoidal and satisfy the boundary 
conditions. Likewise, u s is assumed to be solenoidal, to satisfy the boundary conditions, 
and to be orthogonal to the bases Uj. We assume also that each of the basis functions has 
unit norm, i.e. ||iij|| = 1. For brevity, we denote by S the set of all possible vector fields 
that are solenoidal, satisfy the boundary conditions and are orthogonal to all Uj, so that 
u s G S. 

In order to address the global stability of the nonlinear Navier-Stokes system Q, we will 
partition its dynamics into the interaction of an ODE, representing the evolution of the 
basis weights a, and a PDE, representing the remaining unmodeled modes of the system u s . 
We work initially with the ODE part only, and hence assume initially that u s = 0. 

First substitute ^ into ( 4a ) and take an inner product of both sides with each of the basis 
functions Uj in turn, yielding an ODE in the form 

dj + (uj, Uj • Vu fe ) cijCik + (Uj, S(uj, u)) dj = — (uf, V 2 Uj) ctj. (10) 

Defining matrices A, W and {Q J } . gl n such that 

Ay := (uj, V 2 uj) , Wij := - (ui,S(\lj,u)) , Q\ k := - (u^ uj ■ Vu k ) , 

with L := -^A + W, and defining a linear matrix- valued operator N : W 1 — > M. nxn as 

N(a) := aj Q j , 



one arrives at a compact representation of the ODE ( 10 ) 



a = f(a) := La + N(a)a. (11) 



S 



Two useful general observations about this system are that the matrix A is symmetric 
and negative-definite (in particular, it is diagonal if the basis functions are chosen as 
eigenfunctions of ([7]) with e = 0), and that a T N(a)a = for all a. The latter assertion is 
a restatement of the energy conservation relation (13J) in finite dimensions. 



3. Stability of Finite Dimensional Models using SOS 

For simplicity of exposition, we will assume in this section that the steady solution u is 
spanned by the basis functions Uj, i.e. that there exist some real constants q such that 



u = UjCj. In this case, one can rewrite the dynamics of the finite dimensional system (11) 
in the equivalent form 

a = —Aa + N(a + c)(a + c)-N(c)c, (12) 
Re 

where we have used the identity Wa = N(a)c + N(c)a. Note that a = is an equilibrium 



solution to ( 12 ) , and we wish to find the largest value of Re for which this system is globally 



asymptotically stable. 

To this end, we first recall that an ODE system a = f(a) is stable if one can find a con- 
tinuously differentiable Lyapunov function V satisfying each of the following conditions [61 
Thm 4.1]: 

V(0) = (LI) 
V{a) > Va / (L2) 
V a V{a) ■ /(a) < Va / 0. (L3) 

There is unfortunately no known method to construct such a function for an arbitrary 
system of nonlinear ODEs. However, in the case of a system described exclusively by 



polynomial functions such as (11), the situation is more hopeful. 



First define the energy-like functions Eg : M n — > M as 

E e {a) := l 2 \\a + ec\\ 2 . 

Of special interest will be the perturbation energy function Eq and the total energy func- 
tion E\. In particular, a useful observation is that 

V Q £ (a) • N(a)a = a T N(a)a = 0, 



i.e. the nonlinear part of the dynamics of the system ( |11[ ) is invariant with respect to the 
perturbation energy. Selecting as a candidate Lyapunov function V = Eq, stability of the 



system (11) is therefore assured for all Re such that 



V a V{a) ■ f = a T La = a J (\-K + W)a < Vx / 0. (13) 

R 
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Calculation of the maximum value of Re for which (13) holds is then straightforward, since 



one needs only to find the largest Re such that the matrix 2A + Re ■ (W + W T ) remains 
negative-definite. Of course, this mirrors exactly the situation in the infinite dimensional 



case. 



We next consider whether it is possible to establish stability of the system (11) using some 



alternative polynomial Lyapunov function. In order to restrict the overall size of our search 
space, we first consider the essential features of such a function. 



3.1. System behavior for extreme values of \\a\ 



Consider first the linear part of the system (11) in isolation, i.e 



a = La. 



(14) 



If L has any positive eigenvalues then the system (14) is unstable, implying immediately 



that the nonlinear system (11) is also unstable. If the system (14) is asymptotically stable 



then there exists some P y such that V(a) = aJPa is positive-definite and 

V a V(a) • f(a) = a T (L T P + PL) a < Va / 0, (15) 
see [6l Thm 4.6]. Such a function also ensures stability of the nonlinear system (|11[) for 



some region around the origin, since the linear component of ( 11 ) dominates when ||a|| <C 1 



Considering the nonlinear term of (11) in isolation, one typically expects that for any 
V(a) = a T Pa with P y 0, 

{a | V a y(a)-iV(a)a>O}/0, 



unless P oc I. Since the nonlinear component of ( |11[ ) is the dominant term when [|a|| 3> 1, 
we should generally not expect to find a second-order positive-definite polynomial Lyapunov 
function V other than the perturbation energy function V = Eq, or some monotone function 



thereof. On the other hand, using the system representation (12) it follows that 

1 



V«#i(a)./(a) = (a + c 
1 



Re 



Aa + [N(a + c)} (a + c) - N(c)c 



Re 



(a T Aa + a T Ac) - a T N(c)c. 



(16) 



Consequently, V ' a E\ • / < for all ||a|| sufficiently large with respect to a fixed Reynolds 
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numbei|^J though the choice V = E\ would not satisfy condition (LI). A reasonable 
approach therefore is to search for a candidate Lyapunov function in the form V = A + B, 
where the components A : M. n — > M. and B : W 1 — > M. have the following properties: 







[A + B](0) 

V a [A(a) + B(a)]^a T P V||o||<l 

V a B(a) « a(a) ■ V £d(a) V||a||>l 

||V a A(a)|| < ||V«B(a)|| V||a||»l, 



(17a) 
(17b) 
(17c) 
(17d) 



where P y and a : W 1 — > R is a nondecreasing positive function. The condition (17b) 
ensures that V a y • / satisfies approximately the linear Lyapunov condition (15) in a local- 



ized region about the origin. The conditions (17c)-(17d) ensure that V a V ■ f < for all 



states sufficiently far from the origin, in accordance with ( 16 ) 



In order to exploit SOS techniques, we restrict our attention to cases where both A and B 
are polynomial functions and degA < degB. A useful observation is that any choice of B 
in the form 

u 

1 



k 

B(a) = T[E e .(a) with J = 1 and 9 

i=l 







(18) 



with k S N satisfies the condition (17c). In searching for a Lyapunov function in the form 
V = A + B, we will view the function A as a term to be optimized, and therefore refer to it 
as the variable term. We will restrict the function B to be some combination of energy-like 



functions in the form (18), and hence refer to it as the energy term. 



3.2. Lyapunov Function Generation Using Sum- of- Squares 

If we restrict our attention to polynomial functions V with no constant term (so that 



V(0) = 0), then the Lyapunov conditions ( |L1| )-(L3) can be rewritten as 

{o | V(a) < 0, l x {a) + 0} = 

{a | -V a V(a) ■ f(a) < 0, f 2 (o) / 0} = 0, 



(19a) 
(19b) 



This effect is not exclusive to the total energy function E\. If one defines the energy-like function 
Ed := |(a + d) T (a + d), then the quadratic part of X7 a Ed ■ f is negative-definite whenever the Reynolds 
number Re and vector d are contained in the set 



(Re, d) 



A + 



Re 



W + W' 



W(d) + W J {d)^ -< 



where W(d) is linear in d and defined such that W(d)a = N(d)a + N(a)d. The above set is convex in Re 
for fixed d and vice- versa. In the case that u = CjUj, in follows that W = W(c) and one can make the 
particularly convenient choice d = c, so that the above set is unbounded in Re. 
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where positive-definite polynomial functions l{ are used in place of the vector-valued con- 
dition a ^ 0. For simplicity, we can define the functions 1% as 



and impose a strict positivity constraint on the values e^. Straightforward application of 
the Positivstellensatz (see [13] and the references therein) shows that satisfaction of the 



conditions (19) is assured if one can identify polynomial functions (s\, s 2 ) such that 



V(a) - lx{a) = si, si G S n 
-V a V(a)-f(a)-£ 2 (a) = s 2 , s 2 eS n , 

where E n denotes the set of all sum-of-squares polynomials in M n 



(SOS) 



The problem of determining whether (SOS) can be satisfied can be reformulated as a 



convex optimization problem in the form of a semidefinite program (SDP) using standard 
software tools El El]- If deg V = 2d (note that the degree of V must be even for (LI ) to 
be satisfied), then the general form of our problem is: 



(SDP) < 



min 

K,Hi,H 2 ,{ey} 

subject to: 



dV 
da 



V(x) - h(a) = m d (a) Y Him d (a) 
■ f(a) - h{a) = m d (a) 1 H 2 m d {a) 

(ei,j,e 2 j) > e Vj G {1, . . . , n}, 



(20a) 

(20b) 
(20c) 

(20d) 
(20e) 



where m d {a) is a vector of all monomials in a with degree less than or equal to d. The 
objective function in our optimization problem is zero since we are interested only in feasi- 
bility. Note that any solution to the problem (SDP) will satisfy the original sum-of-squares 



condition (SOS), since the semidefiniteness constraint (20d) ensures that (20b)-(20c) can 



be expressed as sums-of-squares following a suitable similarity transformation. The lower 
bounding constant e for the terms e{j in (20e) must be strictly positive, though it is other- 
wise arbitrary. 



3.3. Determining Stable Values for Re 



One can estimate an upper bound on the value of Re for which a solution to (SOS) can 
be found via a straightforward binary search strategy. However, in all but the trivial case 
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V = Eq, there is no reason to suppose a priori that if a solution to (SOS) can be found for 
some Re, then a solution can be found for all Re £ (0 Re]. Provision of such an assurance 
is possible by augmenting (SOS) with additional constraints. First note that the Lyapunov 
condition (L3) can be written as 



V a V(a) ■ Wa + N{a)a + 



Re 



V a V{a) ■ Aa < Va / 0. 



(21) 



If (21) is satisfied for some Re, then it is satisfied for all Re G [0, Re] provided that the 
second term V a V(a) • Aa < for all a. Satisfaction of this condition can be imposed as a 
sum-of-squares constraint, 



V a V(a) ■ Aa = s 3 , s 3 G S n , 



(22) 



and included as an additional condition to (SOS) (or checked a posteriori). 

Given a Lyapunov function V for some value Re, it is possible to compute directly the 
smallest and largest value Re for which V is a Lyapunov function, since (21) is affine in 
1/Re; e.g. one can compute an upper bound by solving the sum-of-squares problem 

I /Re 

1 



mm 

i?e>0 



subject to: 

and taking the inverse of its minimum value 



V a V{a) (Wa + N(a)a) - — (V a V(a) ■ Aa) - £ 2 (a) G S w 



3.4- Computational Complexity 

We next consider the computational effort required to solve the problem (SDP) for various 
degrees of Lyapunov candidate function V. If we assume that V is a polynomial function 
with arbitrary coefficients and deg V = 2d, then the monomial vector m(a) is composed of 
D distinct monomial terms, where 

(n + rf)! 
raid! ' 

Standard results from semidefinite programming ensure that one can solve the problem 
(SDP) in 0(t/~D) iterations using a primal-dual interior point methocQ with each itera- 
tion requiring O(D^) operations. In practice, it is generally the case that the number of 



4 More precisely, one can guarantee that a primal-dual interior point algorithm will reduce the duality 
gap of its solution iterate to a multiple e of its original value within 0(ln(l/e)\/73) iterations. The reader 
is referred to [2[ 1181 119] and references therein for an overview of algorithms and complexity results for 
semidefinite programming. 
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iterations required to solve a semidefinite programming problem is roughly constant with 
respect to increasing problem size, so the computation time is determined almost entirely 
by the per-iteration computation cost. 

The rapid increase in computational burden with increasing system dimension means that 
SOS methods are likely to be applicable for relatively low dimensional models only, even if 



one assumes that the considerable degree of problem-specific structure inherent in ( 20 ) can 



be somehow exploited (e.g. using a structured approach such as (17)). In particular, it is 
not advisable to attempt to estimate the maximum stable Reynolds number in the infinite 
dimensional Navier-Stokes system Q via solution of a succession of problems in the form 



(20) with increasing dimension. We therefore require a more indirect approach, whereby 
the finite-dimensional techniques of this section can be extended to the infinite-dimensional 
system Q without excessive additional computation. We propose such an approach in the 
remainder of the paper. 



4. Infinite Dimensional Flow Models 



We now return to the general case where u s ^ 0, which we will view as an uncertain forcing 
term in our ODE. In this case substituting ^ into (4a) and taking an inner product of 
both sides results in a model similar to the ODE (11), but with additional perturbation 
terms in u s , i.e. 

a = f(a) + e a (u s ) + Q b (u s , a) + c (u s ), (23) 
where the additional perturbation terms are defined as 

@ ai (u s ) := — (iij, V 2 u s ) - (ui,5(u s ,u)) (24a) 

0&i(u s , a) := - (ui,S(uj, u s )) aj (24b) 
© ci (u s ) := - (iij, u s • Vu s ) , (24c) 



and /(a) is as defined in ( pTT[ ) . In the above, a subscript i indicates that the expression 
is the ith element of a vector quantity. The perturbation term O a represents a linear 
disturbance in u s , 6;, represents a bilinear disturbance in (u s , a) , and O c represents a 
quadratic disturbance in u s . 

We would like to bound the influence of each of these perturbation terms on our ODE 
in terms of ||u s || and [|o||. In order to do so, we apply ([2]) repeatedly to eliminate the 
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appearance of terms Vu s , so that (24) can be rewritten as^] 

©ai(Us) = (u s ,hj) , hi := — 

Re 

0&i(u s , a) = (u s , hij) dj, hij := u, ■ Vui - Vu^ • Uj, 
@ci(u s ) = (u s ,u s • VUj) . 



V 2 Uj — Vu • Uj + u • Vui, 



(25a) 

(25b) 
(25c) 



We are of course left with an ODE in the form ( 23 ) which still features the perturbations u s . 



We next bound the influence of this term by modeling only the evolution of its energy q, 



which we model as q = ||u s || /2. In the process we add a single ODE to supplement (23) 
representing the time evolution of the squared energy term q 2 . 



Substituting ^ into (4a) and taking an inner product of both sides with the total velocity 
field u = Ujaj + u s provides the additional ODE in term of the perturbation energy q 2 , 



where 



(<7 2 ) = o T /(o) - a J a + r(u s ) + x(u s , a), 



r(u s ) := — <u s , V 2 u s ) - (u s , 5(u s , u)) , 



X(u s ,a) := (u s 



Si 



Re 



(26) 

(27) 
(28) 



Verification of the above relies on the aforementioned assumptions about the subspace 5 
and on application of the various identities described in Section [T] In particular, these 
allow one to establish the relations 



u 



u 



d_ 
' M 



a 1 a + (q 2 ) and a T /(c 



1 

Re 



[ui, V 2 uj) - (ui,5(uj,u)) 



Note that in (26), the terms a T f(a) and T(u s ) represent the self-contained dissipation or 
generation of energy depending on a^Uj and u s , while the term x(u s ,a) represents the 
generation or dissipation of energy containing cross terms between these velocity fields. 



4-1. Description as an Uncertain System 

The complete system of interest can now be written as 

a = f(a) + a (u s ) + Q b (u s ,a) + O c (u s 
(g2) = a T /(a) - a T a + r(u s ) + x(u s , a). 



(29a) 
(29b) 



'The notation used can be clarified by the equivalent expression for the Cartesian components of the 



vector h»: ftf = ^VV 



Ou „ k i — V7, T) 

q^Ui + u • VUi 
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We are now free to treat u s as an uncertain term driving the ODE system (29), whose 
time evolution is known to satisfy the subspace constraint u s £ S and the energy con- 
straint q 2 = ||u s || 2 /2. The worst-case effect of this uncertainty can then be bounded via 
appropriate norm bounds. 



The first of these bounds relates to the uncertain terms in (29a). There exist constants 
d > and a polynomial function pi(a, i 



> such that 

2 



|6 a (u s ) + Q b (u s ,a) + 6 c (u s )|| < pi(a,q) 



ciq 2 + c 2 q 2 \\a 



| 2 W 



(30) 



for any a and u s . A rigorous proof of the existence of these constants is given in Appendix 



|A| Critically, estimation of the coefficients q involves the solution only of linear problems 
for partial differential equations and optimization over finite-dimensional polynomials. 

A second bound relates to the uncertain term r(u s ) in (29b). Comparing (27) with ^ 
shows that r(u s ) < k ||u s || 2 with k = X±. However, u s satisfies the additional constraints 
(u s , Uj) = and therefore may admit a stronger bound. Note that the number of positive 
eigenvalues of ([7]) is always finite [I]. Hence, if are chosen as the first n eigenfunctions 
of ([7]) and n is large enough, then 



rfu.l < k.c 



U. s 



2n s q 



for all u s G S, where k s = \ n +i < 0. If Uj are not eigenfunctions of ([7]), then k ( 
largest eigenvalue of the following problem 



(31) 
is the 



Au + /i fc u fc = (e 
(u fc ,u) = 0, V 



1 



Re 

u = 0, u\ m 



A)u + X7p 
0. 



In what follows we will assume that k s < in (31). 



A final bound relates to the uncertain term x( u s) i n (29b). If Uj are eigenfunctions of 
([7]) then x = 0, because in this case gj = — AjUj + V0i with some scalar functions <f>i 
and because u s is orthogonal to both Uj (by definition) and to gradients of any scalars 
(since V • u s = 0). In the general case there exists a constant d and a polynomial function 
P2(a, q) > such that 



|x(u s ,a)|| <P2(a,q) 



7 2 llall 2 . 



(32) 



The proof is very similar to the proof of (30). 



5. Stability of Infinite Dimensional Models using SOS 



Given the (uncertain) ODE system (29), we can now search for a Lyapunov function 
verifying stability of the composite state vector (a, q 2 ). We therefore would like to construct 
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a Lyapunov function V : M n x R — > E such that 

+ W) ' (g2) < °' v(a ' Us) Us e 5 ' (33) 

We can expand the left hand side of this condition and collect terms to get the equivalent 
Lyapunov condition 



dV 
da 



. dv _ fdv dv T \ /_ _ \ ay 



(34) 



where we have omitted the arguments for (/, T, a , 0&, C , x) for brevity. For simplicity, 
we will assume that the function V is chosen in such a way that 

dV/d{q 2 ) > (35) 

Consequently, 



dV 



, r(u s ) < 0, Vu s G 5\{0} 



provided that (31) and (35) hold. 



5.1. Comparison to the Energy Method 

Note that if one chooses a candidate Lyapunov function by making the most obvious 
generalization of the type of function suggested in Section [3J i.e. if one chooses V(a, q 2 ^ 



V a (a) + Yli = i(Ee i (a) + q ), where V a {-) is some polynomial function, then (35) is satisfied. 
The term 

fdV dV n 
\ da d(q 2 ) 



in the Lyapunov condition (34) can be viewed as a misalignment between the (scaled) 
gradient of the energy function Eq{o) and the gradient term dV /da. If one chooses as a 
candidate Lyapunov function 

V(a,q 2 ) =E (a) + q 2 , 

then the above misalignment term is zero. If, additionally, one chooses Uj such that 
X = 0, which was shown above always to be possible, the situation reduces to the usual 
global stability condition using energy functions. Consequently, if energy can be used as 
a Lyapunov function for the system Q for some Reynolds number Re, then the choice 



V(a, q ) = Eq{o) + q will satisfy the conditions (34). 



When the system remains globally stable for Reynolds numbers beyond this energy stability 
limit, one should first ask whether there exists any polynomial in (a, q 2 ) that will serve as 
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a Lyapunov function. We give a constructive proof of the existence of such a function in 



Appendix B 



It remains to demonstrate that the Lyapunov function V satisfying ( 33 ) can be constructed 
in a systematic way using the SOS approach. 



5.2. Conversion to a Sum- of- Squares problem 



After applying (31), the inequality in the Lyapunov condition (34) can be written in vec- 
torized form as 



-d-J + W)' 2Ksq]< 



dV dV _T\ dV 

" 1 ' st<F) 



da djq 1 ) 1 



X 



(36) 



We next apply the Schwarz inequality, ( 30 ) , ( 32 ) , and ( 35 ) to arrive at a sufficient condition 



for satisfaction of the inequality ( 33 ) : 
dV 



dV 

da f{a) + W)' 2Ksq2)< ~ 



dV dV „T\ dV 
da dj?) a I ' djq 1 ) 



[Pi (a, q) + Vi (a, q)] 2 , V(a, g) / 0, 

(37) 



where | • | represents the standard Euclidian norm in IR n . The above can be rewritten more 
compactly as 

g(a,q)<-\\h(a,q)\\-ph(a,q), V(a,g)^0, (38) 



where 



, , dV 3V 2 

g(a,q) := —Ka) + ^.2K.q 

p(a,q) :=px(a,q) +p 2 (a,q). 



We next apply the following matrix property, based on the Schur complement [21 Sec. A. 5. 5] . 
For any vector u and scalar t, 



u < t 



t u> 
u tl 



>- 0. 



The condition (36) is therefore equivalent to 



g(a,q) h T (a,q) ■ p*(a,q) 

h(a,q) ■ p2(a,q) g(a,q)I 



-<0, V(a,g)/0, 
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Since p is a nonzero function for all nonzero (a, q), the above can be pre- and post-multiplied 



by ( o ) to get the equivalent matrix inequality condition 



H(a,q) :-- 



g{a, q) ■ p{a, q) h T (a, q) ■ p{a, q) 
h(a,q) -p(a,q) g(a,q)I 



-<0, V(a,g)^0. 



(39) 



The most important thing to note about (39) is that it is linear in the coefficients of 



the Lyapunov function V. This linear matrix inequality (LMI) can be converted to an 
equivalent scalar polynomial inequality via introduction of an additional variable z G M n+1 . 
It is straightforward to verify that (|39|) is equivalent to 

(40) 



z T H(a, q)z < 0, Vz / 0, V(a, q) / 0. 



In other words, if the function V is chosen such that dV/d(q 2 ) is nonnegative and 



[a,q,z) 



and 



z T H(a,q)z>0 
\a\\ 2 +q 2 ^ 0, \\zf / 



(41a) 



(a,q) 



V(a,q 2 ) < ' 
\\af + q 2 ^0 



(41b) 



then the Lyapunov condition (34) is satisfied. This is a standard form convenient for 



applying the Positivstellensatz theorem [14] . from which it follows that (41) are satisfied 



if and only if there exist non- negative integer values Mi , M2 , M3 , M4 and M5 and sum-of- 
squares of polynomials Sj(a, q, z) and <Tj(a, q) such that 



Afi 



so 



Y J ^[z T H(a,q)z)Y-(\\a\\ 2 + q 2 ) 
i=l 



2\2M 2 || Z || 4M 3 



and 



M 4 



|a|| 2 + (? 2 ) 2M5 . 



vo = -^2aj[V(a,q 2 W 
J'=l 

One can now choose (from empirical considerations or by trial and error) the integers Mj 
and the polynomials Sj and <jj for > 0, thus obtaining expressions for sq and <7o via 
the coefficients of the polynomial V. Then determining the coefficients of V such that so 
and do are sum of squares can be attempted using the existing packages SOSTOOLS [12] or 
YALMIP 0[B]. For the purposes of computational efficiency, the selection of V should be 
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subject to structural constraints similar to the case of an ODE system obtained by simple 
truncation, as in Section [3} 

Finally, we note that if one is able to choose the bases Uj such that x = 0, then the 
preceding problem can be simplified somewhat. In this case, one is free to define 

h(a,q) := - T^^ya 1 "^ , p(a, q) := pi(a, q) , 



in (36)-(38). It is then easy to show that the choice of Lyapunov function V = Eq{o) + q 



will satisfy not only the stability condition (34), but also the robust LMI condition (39) 



whenever the total perturbation energy is a Lyapunov functional for the Navier-Stokes 
system Q. This ensures that the proposed method will always yields results at least as 
good as classical energy perturbation methods. 



6. A Finite Dimensional Example 

In this section we present numerical results for a model of Couette floftj^] using the finite- 
dimensional stability analysis results of Section [3| Couette flow refers to the shear flow of 
a fluid between two infinite parallel plates as shown in Figure [T] 




Figure 1: Periodic flow between parallel plates. 

We employ the finite-dimensional ninth-order model developed in [9] for this flow, and 
make assumptions identical to those in [3 [TU] for the purposes of comparison. The volume 
force is assumed to be _ 

f = (^ S in(vry/2);0;0), 

and the flow is assumed to be periodic in the spanwise and streamwise directions, with 
wavelengths L z and L x respectively. We fix L x = 4tt and L z = 2ir throughout, with plate 
separation h = 2. 



The results of this section were first reported in [4]. 
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The flow is assumed to have free-slip boundary conditions 



ly\y=±l 



0, 



du x 
dy 



y=±l 



duz 
dy 



0. 



y=±l 



A set of nine basis functions Uj were selected in [9] based on physical insights and observa- 
tions arising from numerical simulation and experiment. For reference these basis functions 



and their expansion into a nonlinear ODE in the form ( 12 ) (equivalently ( 11 )) are included 



m 



Appendix C . In this example, ui is a laminar solution to the Navier-Stokes equation ([!]). 



We compute an upper bound on the value of Re for which the system (11) is guaranteed 



to be stable using a variety of Lyapunov candidates, each of which satisfies the structural 



conditions (17)-(18). 



To compute the upper bound, we use a bisection method to find the largest Reynolds num- 
ber for which the sum-of-squares optimization problem (20) could be solved with e = 1CP 5 



in (20e). All of the results obtained were computed on a 2.33 GHz Intel Xeon processor 



with 3.6 GB RAM, using the YALMIP interface to the SDP solver SeDuMi EJ. 

Overall results are summarized in Table [TJ In each case we provide the form of Lyapunov 
function used and the maximum value of Re for which a Lyapunov function in this form 
could be identified. We also report the total solver time required (which includes the time 
spent both in the SDP solver SeDuMi and in preprocessing tasks by YALMIP), the total 
number of monomial terms in the vector m^(a) that appears in the equality constraints of 



(20b)-(20c), the number of decision variables in V that take nonzero values in the solution 



to (20), and the total number of nonzero elements required in the solution for the matrices 
Hi and H%. 

The largest Reynolds number for which a Lyapunov function could be identified was Re = 
54.1, which compares very favorably to the value Re = 7.5 for which a purely energy-based 
method succeeds. Previous results from numerical work in [10J, using the same model, 
have suggested a value of no more than Re ~ 80 before the system becomes unstable. 
This suggests that the method we propose is not unduly conservative. Note that all of the 
computed stability bounds on Re presented in Table [T] are unchanged if one includes the 
additional constraint (22). 



6.1. Perturbation Energy as a Lyapunov Function [Case 1] 



We first consider the simplest case where one takes V = Eq. Recalling (13), an upper 



bound on the value of Re for which the system is guaranteed stable is readily found via 
solution of the following semidefinite programming (SDP) problem: 

max Re 

s.t. 2k + Re-(W + W T ) ~< 0. 
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Case 


Lyapunov Function 




solver 
time (sec) 


number of 
monomials 


nti7i Pi ryr 

llllZj 1 J I Ul 

dim(p) 


nri7f ffi \4- 

llil/j y 1 1 1 JH^ 

nnz(H 2 ) 


1 


V = E 


7.5 


0.1 








2 


V = \a T Pa + E E 2 


23.9 


2.7 


54 


21 


1276 


3 


V = §a T Pa + E EiE 2 


28.5 


20.8 


219 


21 


19410 


4 


V = m 2 (a) T Pm 2 (a) + E E 1 E 2 


54.1 


43.3 


219 


776 


24042 


5 


V = p T m 4 (a) + EqExE 2 


54.1 


41.4 


219 


190 


24042 



Table 1: Computed lower bounds on maximum stable Reynolds number for various Lyapunov functions 



This method is analogous to the use of the conventional energy-based approach described 

6.2. Lyapunov Functions with Second- Order Variable Terms [Cases 2 & 3] 

We next consider candidate Lyapunov functions with the variable term A restricted to the 
quadratic form A(a) = a? Pa, where P G S n is taken as a decision variable to be optimized, 
and B restricted to a weighted sum of energy terms of higher order. We restrict P to those 
matrices whose sparsity patterns match that of solutions to the Lyapunov equation 

L T P + PL = -I, (42) 

and note that the sparsity pattern of such solutions are invariant with respect to Re. 
This choice does not result in any apparent increase in conservatism (i.e. it does not affect 



the range of values for which (SOS) can be solved), while reducing considerably the overall 
computation time required. Increasing the degree of the energy term between cases 2 and 3 
shows a slight improvement in the maximum value of Re for which stability can be assured. 

6.3. Lyapunov Functions with Fourth-Order Variable Terms [Cases 4 & 5] 

Finally, we consider candidate Lyapunov functions with variable terms of fourth order. 



In these cases, direct solution of (SOS) requires substantially increased computational 
effort relative to cases with second-order variable terms. For case 4 we take A(a) = 
m,2 (a)" 1 " Prri2{a) with the matrix P treated as a symmetric decision variable. For case 5 
we take A(a) = p T rrii(a) with vector p treated as a decision variable. In both cases, A is 
restricted to contain only terms of degree at least 2. 



In both cases, a solution to (SOS) was found once initially for an arbitrarily chosen (small) 
value of Re in order to identify a likely sparsity pattern for P and p respectively. Subsequent 
computations enforced sparsity of P and p by setting to zero those elements taking relatively 
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small values (i.e. < 10~ 7 ) in the first trial. In both cases, this procedure resulted in a 
substantial reduction in the degrees of freedom afforded to the solver, with a consequently 
large reduction in overall computation time. 

A selection of cross-sectional plots showing the phase space of the system a = /(a) = 
L(a) + N(a)a and level sets of the computed function V for Case 5 are shown in Figure [2] 
The following features are of interest: 



Figures 2|(a) and 2|(c) show clearly that level sets of V are centered on a 
\\a\\ <C 1, and are centered on a = c = (1, 0, • • • , 0) for ||a|| 3> 1. 



for 



Level sets of V are radially symmetric in 2|(b) Streamlines of / approach tangents 
to these level sets as Re — > 00. Radial symmetry is also clear in 2|(d) for ||a|| 3> 1. 



There is significant warping of the level sets of V (with respect to a radially symmetric 



energy function) in 



and |2f(d~) 



Optimizing V over a high-order polynomial allows 
increased freedom to shape these sets, thereby increasing the range of values Re for 
which / can be proven stable. 



6.4- Lyapunov Functions with Higher-Order Variable Terms 



Direct solution of problem (20) for Lyapunov functions with variables terms of degree 
greater than four is more problematic given the long computation times required. How- 
ever, it is possible that close scrutiny of the results from cases 4-5 may give some indication 
of appropriate sparsity structures that may be exploited. For example, as shown in Fig- 
ure [3j the sparsity pattern of the matrix P featuring in the term m2(a) T Pm2(a) in case 
4 can be reordered to block diagonal form. Note in particular that the upper-left hand 
corner of the unordered matrix in Figure [3j which corresponds to the second-order terms in 
m,2(a) T Pm2(a) , adopts an identical sparsity pattern to the solution of (42). The resultant 
reordering is such that one can rewrite the variable component of V as 

T 



m2(a) T Pm2(a) 



ni2 

ml 

m\ 
~ 4 



a) 

[a) 



Pi 



Pa 



Pi 





~m\{a)~ 




m 2 2 (a) 








.m|(a)_ 



where each of the matrices Pi is symmetric and dense and the reordered and partitioned 



monomial terms to 2 (o) are defined as 



7712(a) 

7712(a) 
m|(o) 
?74(a) 



(02, a 3 , aia 2 , aia 3 , a 4 a 6 , a 5 a 6 , a 4 a 7 , a 5 a 7 , a 4 a 8 , a 5 a 8 , a 2 a 9 , a 3 a 9 ) T 

(a 4 , a 5 , a 4 a 4 , aia 5 , a 2 a 6 , a 3 a 6 , a 2 a 7 , a 3 a 7 , a 2 a 8 , a 3 a 8 , a 4 a 9 , a 5 a 9 ) T 

(ae, a 7 , a 8 , a2a 4 , a 3 a 4 , 0205, a 3 as, aia6, aia 7 , a 4 a 8 , a6a 9 , a 7 a 9 , a 8 a 9 ) T 

(1, a\, ag, a2a 3 , a 4 a5, a6a 7 , a6a 8 , a 7 a 8 , aia 9 , a 1 ,a 2 , a 3 , a 4 , a 5 , a 6 , a 7 , a 8 , a 9 ) . 



(43) 
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-3-2-10123 -3-2-10123 

X 1 x 6 

(a) Cut through (01,03). (b) Cut through (06,07). 




-3-2-10123 -3-2-10123 



X 1 x 2 
(c) Cut through (01,09). (d) Cut through (02,09) 

Figure 2: Selection of cross-sectional plots of the phase space for the nonlinear system a = La + N(a)a with 
Re = 54. Each figure shows a cut through the origin of the phase space for some (a,i,aj) pair. Shading 
indicates the relative value of the Lyapunov function V computed for Case 5 in Table [T] with dashed black 
lines indicating level sets. Streamlines of the nonlinear system, projected onto the cutting plane, are shown 
in solid blue. 
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sparse 

y 30 

reordering 



50 :::::::::::::::::: 

10 20 30 40 50 10 20 30 40 50 

nz = 776 nz = 776 

Figure 3: Case 4: Reordering of monomial terms in 1112(0) reveals block-diagonal structure in P. 
7. Conclusions 

A new method for analyzing the global stability of a fluid flows has been proposed. This 
method requires only the solutions of linear eigenvalue problems for systems of linear partial 
differential equations, combined with a nonlinear analysis of a system governed by nonlinear 
ordinary differential equations, which can be treated using the polynomial sum-of-squares 
approach. The method is proven always to yield results that at least as good as classical 
energy methods, in the sense that if the global stability of a particular flow can be proved 
by the energy method then it also can be proved by the new method. 

Moreover, if the flow remains globally stable for values of the Reynolds number in a range 
extending beyond the maximum for which global stability can be proved by the energy 
method, then a polynomial Lyapunov function is still guaranteed to exist in at least part 
of this extended range. The methods proposed in this paper can then apply, provided that 
the appropriate quantities can be expressed as sums-of-squares. 

Application of our method to a finite dimensional model system suggests that using the pro- 
posed method might allow proving global stability for much higher values of the Reynolds 
number than the values for which this can be done by the energy method. It would of 
course be of interest to extend such an example to the infinite dimensional case using the 
techniques of Section [5] 
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Appendix A. Calculating the norm bounding coefficients Cj 



In this appendix we will prove (30). The claim is that there exist constants (ci, 02,03) > 
such that 

||9 a (u s ) + e b {u s ,a) + G c K)f < Cl q 2 + c 2 ||a|| 2 q 2 + c 3 q 4 Va G R n , Vu s G S. (A.l) 
First note that for any norm, 

||u + v|| 2 < ||u + v|| 2 + \\u — v\\ 2 = 2 (||ii|| 2 + |M| 2 



We can apply this inequality to (A.l) and then compute bounds for each of the terms in 
turn. 

Part I - Computing a bound on ||0 a (u s ) + 0{>(u s )|| 2 : 

||G a (u s ) + 6 fe (u s )|| 2 = ^2 (u s ,hj + hijOj) 2 < 1(a) \\u s \\ 2 

i 

where 

ti \ Ei ( u s,hi + hijaj} 2 

1(a) = sup 2 — — • (A-2) 

Us \\ U s\\ 

Let hj = h, + Vi^j + /?jfcUfc be a projection of hi on the solenoidal subspace orthogonal to all 
Uj, so that Vhj = and (hj,Ufc) = 0. Then the scalar fields 0, should satisfy the Poisson 
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equation V 2 <; 
condition 



-V • hj. We will make the solution unique by imposing the boundary 







n 



-(hi,n), 



no 



where n is a unit vector normal to the boundary <9f2. This condition ensures that the 
normal component of hj equals zero: (hj,n) = 0. Similarly, we also introduce solenoidal 



fields hjj-a. 



Strictly speaking, in (A. 2) the solenoidal velocity field u s should also satisfy the full bound- 
ary condition u s \ dn = 0. However, this condition can be relaxed to (u s ,n)| 9n = without 
changing the value of 1(a), since for any vector field v satisfying only (v,n)| 9n = one 
can easily construct a solenoidal field v' satisfying v'\qq = and such that (v — v', hj) 
arbitrarily small. With this relaxation one can represent u s in (A.2) as 



is 



Us 



fjjhj + Uijhij + u s , 



where u s is orthogonal to all hj and . Obviously, u s does not contribute to the numerator 



of (A.2) but can only increase the denominator. Therefore, for calculating the supremum 



one can take u s = 0. Then (A.2) becomes 



1(a) 



sup 

u k ,u kl 



Ej (Ukh-k + U M h kl , hi + hijttj 



Ukhk + Ukih k i 



(A.3) 



Since 



yUkhk + Ukihkl, hj + hy%\ (hi + hijaj, hj + hjjOj 



U k h k + Ukih 



ki 



hj + hijaj 



both the existence of the supremum and of the coefficients c\ and C2 such that 1(a) < 
(c\ + C2 ||a|| 2 )/2 becomes obvious. Note that the numerator of the estimate (A.3) is a 



fourth-order polynomial in coefficients Uk, Um, and aj, and that polynomial is quadratic 
separately in U^, Um, and aj, which implies that efficient numerical methods can be found 
for determining c\ and C2- 



Part II - Computing a bound on ||0 c (u s )||: 

Calculation of this bounds is straightforward, and amounts to finding a solution to a 
generalized eigenvalue problem using a symmetric version of the operator Vuj , in a manner 
similar to that used to calculate the bound A in ([7]). This results in a constant c%i such 
that 



|Oci(u s 



< 



cm 



so that C3 = 2 c ii) provides a conservative bound. This completes the proof of (A.l). 



28 



Appendix B. Existence of a polynomial Lyapunov functional for Reynolds 
numbers greater than the energy stability limit 



In general a flow can, and often will, remain globally stable in a certain range of Reynolds 
numbers greater than the maximum Reynolds number Re e for which global stability can 
be proved by the energy method. We will prove now that at least in some part of this range 
there exists a Lyapunov functional that is polynomial in a and q, i.e. of a form suitable 
for the proposed method. To this end we will simply give an explicit expression for a 
function which is polynomial in a and q 2 and which satisfies all the conditions for Lyapunov 
functionals. For this purpose we will assume that Uj are chosen as the eigenfunctions 
of ([7]), and note that this is always possible. As a result, we will have x = as explained 
in Subsection |4.1| We will, of course, use explicitly also the assumption that the flow is 
indeed globally stable in some vicinity of the energy stability limit. More precisely, we 
require that 

n 

£/?/[(l + 2£ )£o] (B-l) 

i=2 

is strictly positive when ai = as = ■ ■ ■ = a n = for all values of a\ including zero and 
infinity. If this were not true at some finite 01 = a± s for any n that would mean that ai s ui 
is a steady non-zero solution, which contradicts the assumption of the basic flow being 
globally stabltQ This means that this condition can be satisfied by selecting large enough 



n. To ensure the positiveness of (B.l) at the origin we will additionally require that the 



flow is linearly asymptotically stable, which is not much of a loss of generality for the flow 



which is already assumed to be globally stable. To ensure the positiveness of (B.l) with 
a\ at infinity, we will also assume that (uj, ui • Vux) ^ for at least one % in 2, . . . , n. If 
ui • Vui 7^ 0, which is often the case, this assumption can always be ensured to be true by 
selecting a sufficiently large n. Finally, we will also assume that the eigenfunctions Uj and 
eigenvalues Aj are continuous functions of Re in some vicinity of the energy stability limit. 

The basic idea is simple. At Re = Re e energy satisfies all the conditions for being a 
Lyapunov function everywhere except along m, where its time derivative is zero. We add 
to it a small correction constructed in such a way that its time derivative is negative along 
ux- The modified expression will therefore be a Lyapunov function at Re = Re e and, by 
continuity, also in its vicinity. In the sequel, we exploit this idea in a more formal proof. 



Note that this argument is valid only if n > 1, and indeed one can check that for n — 1 our approach 
will not work for Re > Re a because the one-dimensional system turns out to be unstable. It can be shown 
that in the general case n has to be greater than the number of positive eigenvalues \i . 
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We propose the following expression for the Lyapunov function 



V = E + q 2 + (E + 



J2\2 



' ^2 aifi 



(B.2) 



i=2 



where f(a) is defined by (11) and e is a small positive value. Note that the sum in (B.2) 
starts at i = 2. Note also that are of order a for small a and of order a 2 for large o. 

We first provide conditions under which this choice of V is positive-definite. As required, 
V = if both a = and q = 0. Otherwise, V > provided e > and 



e < mm 



Ep + -Bp 



This ratio is easily shown to be bounded on an open ball of sufficiently small (large) radius. 
The fraction is continuous outside (inside) of such a ball. Existence of a positive minimizer 
is therefore ensured by the extreme value theorem. 

We next develop conditions under which V is negative definite. Following the same argu- 



ments as in Subsection 5.2 but noticing that here x = by assumption, we arrive at the 



following sufficient condition for satisfaction of the Lyapunov condition ( 33 ) : 
dV 



da 



f(a) 



dV 

W) 



2K s q' 



< 



dV 
da 



dV 

W) { 



Pi(a,q), V(o,g) ^ 0. 



(B-3) 



From the definitions of Eq and \ and from the energy equation ^ it follows that 

i=l 



da 



iof. 



Using also the definition of pi in (30), the inequality (B.3) can be rewritten as 
Xiaj + 2n s q 2 + Y?i=2 X i a t 



E + 



- + 



e £?=2 iff + at ELi +e\§- a *itt IffK^i + c 2 ||a|| 2 + c 3 g 2 )§ 



:i + 2^ + 2g 2 )(B;o + g 2 



< 0. 



30 



If we define 



fi 
Fi 



ai/(E + q z )z, q:=q/(E + q 
f t /[(l + 2E + 2q 2 )(E +q 2 )}K 



2\ i 
Z ^ 2 



Fi(a,q) :-- 



Z^j=l da, h 



(l + 2E + 2q 2 )(E + q 2 )2 



F q :=F q {a,q) 



d_ 
da 



i=2 



(ci + c 2 ||a|| +c 3 g 2 )a 
(l + 2E + 2q 2 )(E Q + q2) l 2 



(B.4) 
(B.5) 

(B.6) 
(B.7) 



then (B.4) can be written more compactly as 

n n 

Aio? + 2n s q 2 + \a\ - e ^ (j? + 5^) + e\q\F q < 0. 



i=2 



i=2 



(B.8) 



Using the identity xy = (Ax) 2 — (Ax — y/(2A)) 2 + y 2 /(4A 2 ), one can rewrite (B.8) as 

n n / n \ 

X 1 ~a 2 + (2k s + eA 2 )q 2 + J> + eA 2 )~a 2 - e £ f 2 + * [f 2 + £ F 2 - 



i=2 



i=2 



i=2 



i=2 / 

S) 2<a < B ' 9 > 



The last two terms of (B.9) are non-positive, so a sufficient condition to satisfy (B.9) is 

n n 

D(Re) = \ x a\ + (2k s + eA 2 + e)q 2 + J^(Aj + eA 2 + e)a 2 - e q 2 + ^(f 2 + a 2 ) 

L i=2 

AA 2 \ r g + 2^A 



i=2 



i=2 



< 0. (B.10) 



Note that if Re = Re e , then Ai = 0, k s < and Aj < 0, with both inequalities strict. 



Satisfaction of (B.10) is therefore ensured if both 

e < -mm{2K s ,\i}/(A 2 + 1) 
and A is chosen sufficiently large to guarantee that 



e+±Cf!+ai)~ 2 U+±F 2 )>* 

i=2 \ i=2 / 
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for all (a, g) ^ 0. A suitable choice is 



A 



F 2 



iC/j' + a?) 



(B.ll) 



provided that this fraction can be shown to be bounded for all (a, 



We first consider whether the numerator of (B.ll) is bounded above. The functions Fj 



and F q are easily shown to be bounded inside (outside) an open ball of sufficiently small 
(large) radius, and are continuous elsewhere. Boundedness over all (a, q) then follows from 
the boundedness theorem. 



The denominator in (B.ll ) is bounded below by a strictly positive value, because we already 
assumed that (B.l) is bounded below by a strictly positive value when 02 = 03 = • • • = 



a n = 0. 

The above argument ensures that D(Re c ) < 0. Then from our continuity assumptions it 
follows that D(Re) < in at least some vicinity where Re > Re e , thus completing the 
proof. 



Appendix C. System Dynamics for Numerical Example in §[6] 

The shear flow model used in the example in Section [6] is taken directly from [9l [10] . We 
include here the basis functions u» and resulting ODE system from [HI [TO] for easy reference. 

The model uses the following 9-dimensional basis of mutually orthogonal, solenoidal basis 
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functions: 



u 3 



U 5 



U6 



U7 



U 8 



7T^ 



'V2sin(vry/2)' 





27 cos(7ry/2) cos(72;) 
7r sin(-7ry/2) sin(7z) 




2 sm(ax) sin(-7ry/2) 

—7 cos(ax) cos 2 (ny / 2) sin(7z) 


asin(ax) cos 2 (ny/2) cos(7z) 

7sin(ax) sin(7ry/2) sin(7z) 


a cos(ax) sin(-7ry/2) cos^z) 

7ra sin(ax) sin(7ry/2) sin(7z) 
2(a 2 + 7 2 ) cos(ax) cos(-7ry/2) sin(7z) 
— 7T7 cos(ax) sin(7ry/2) cos(7z) 



u 2 



114 



cos 2 (-7ry/2) cos(7z) 






cos(ax) cos 2 (7ry/2) 



4 

4 



4^/2 



^3(a 2 + 7 2 )' 
V(« 2 + 7 2 )' 

#8, U9 



V2sm(3vry/2)' 





where a = 2ir/L x , (3 = tt/2, 7 = 2ir/L z and 



iVs 



2^/2 



v / (a 2 + 7 2 )(4a 2 + 4 7 2 + 7r 2 ) 



It is easily verified that w = ui is a laminar solution of the Navier-Stokes equation ([I]) 
when the volume force is 

f= (w sin( ^ /2);0;0 )- 

These basis functions can then be expanded via Galerkin projection into a nonlinear system 
of ODEs as described in Section [2} Define the following notation for neatness: 



K a B := \/a 2 + /3 2 , K ai := \J a 2 + 7 2 , k« 7 := y 7 /? 2 +7 2 , and k q « 7 := -\/a 2 + P 2 + 



Then the ODE in the numerical example is the form (12), with c = diag(l, 0, . . . ,0), 



I 2 4/3 2 1 2 2 3a 2 +4/3 2 2 3« 2 + 4/3 2 + 3 7 2 
A = -diag I p , — +7 , K /37 , , k q/3 , 



, k 2 ^ 7 , k 2 / 3 7 , 9/3 2 
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and 



[N(a 

[N(a 

[N(a 

[N(a 

[N(a 

[N(a 

Ma 

Ma 

[N(a 

where 
arguments. 



a]i 

a] 2 

a] 3 

a] 4 

a] 5 

a] 6 
a] 7 

a] 8 
a] 9 



3 h 



2 K,9 7 

10 7 2 
3\/6 



a2a 3 



3 g7 

2 



-a 6 a 8 



-a4d6 



/fi " a5 ° 7 /fi 

(a 5 a 6 + a 4 a 7 ) + 



3 l^cx'y AC /5"y 



a 5 a 8 



2 a/?7 _ x , ^ 2 (3a 2 + 7 2 )- 37 2 4 7 



a 10 a 2 

^(aia 5 + a 5 a 9 ) - — -= a 2 a 6 

V6 3v6 K «7 



3 a/?7 



/3^7 

2 K ( g 7 

a4a 8 



(aia 3 + a 3 ag) 



76 



(aitt4 + a 4 a 9 ) + 



2 



a 3 a 7 



3\/6 K a7 

a/?7 



-a 2 a 4 



2 a/?7 a' 
a 3 a 6 + -= a 2 a 7 

2 2a/?7 a , 
o 3 a 5 + — =(oia 7 + a 7 a 9 ) + 

O ^07^/37 V 6 



a 2 /? 2 

^ AC a ^AC^^yAC a ^^y 

a/37 



a 3 a 8 



\Z^AC a >yAC a ^ 

3 /3 7 



a 2 a 8 



2 K 



a/37 



(aia 8 + a 8 a 9 ) 



-a 3 a 4 



(-« 2 +7 2 ) 



V^6ACq,^ AC^j-y 

7 2 (3a 2 -/3 2 + 37 2 ) 



-a 2 a 5 7=(aia 6 + a 6 a 9 ) 

VDK Q7 vd 



a 3 a4 + 



2 a/?7 

3 Av Q7 Av a ^3 7 



a 2 a 5 



3^7 

K,3 7 



-a 2 a 3 



3 h 



2 K 



a/37 



-a@a 8 , 



iV(a)a]j is the i th component of N(a)a and diag(-) forms a diagonal matrix from its 
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